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Abstract. We present in this paper a very adapted finite volume numerical scheme for 
transport type-equation. The scheme is an hybrid one combining an anti-dissipative method 
with down-winding approach for the flux [5, 4] and an high accurate method as the WEN05 
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1. Introduction 

In this paper we are interested in finite volume numerical simulations of PDEs of hyperbolic 
type with transport term. More precisely we are looking to correct the numerical diffusion 
which appears in the simulations of the asymptotic profile of such PDEs. Indeed this nu¬ 
merical diffusion is an artifact that is inherent to most of existing numerical schemes even 
for high-order ones. So, as reported first in [1] and confirmed in [5] for the Lifshitz-Slyozov 
equation which is of transport type, capturing numerically the exact asymptotic profile for 
transport equation is a real challenge because numerical diffusion smooths out the fronts 
and leads to an artificial profile. In the context of biology modelling, specially in population 
dynamics, most of existing models [15, 7, 2] contain at least a transport part which takes into 
account the growth of the considered population, it is crucial to recover the exact asymptotic 
profile in order to predict the behavior of the population or to estimate some parameters for 
instance its growth, division or death rates by using measures that are based on the profile. 
Some authors address the question consisting to correct this inherent numerical diffusion by 
establishing adequate schemes such as the WENO (Weighted Essentially Non-Oscillatory) 
scheme [8], anti-dissipative scheme the ADM (Anti Dissipative Method) [5], etc. All these 
schemes define their numerical fluxes in order to minimize at best the artifacts. 

We propose here a hybrid finite volume scheme where the construction of the numerical 
flux appears as a combination of the WENO and ADM fluxes. This scheme is very suitable 
for firstly removing de numerical artifacts but also correct the stair treads appearing in 
ADM scheme. More precisely we use a discontinuity detector at each grid point and use 
the WENO order 5 scheme when the solution is regular near the point otherwise we apply 
the ADM. As a validation of this hybrid scheme we use two kind of test cases. The first 
one is a classical (academic) test case for transport equation where the initial distribution is 
considered in one hand as a very oscillatory one as given in [8]; in the other hand, we use 
in 2D the famous Zalesak’ disk test that is given in [9]. The second kind of test is based 
on population dynamics and polymerization process where we consider a population of cells 
or polymers growing either by nutrients uptake (for cells) or by gain and lost of monomers 
(for polymers). This application on population dynamics is of great importance because in 
many cases some predictions on the numerical behavior of the solution allow to investigate 
inverse problem of estimating relevant parameters of the considered model. So, having a bad 
numerical reconstruction induces bad parameters estimation. 

The paper is organized as follows. In section 2, we recall the biology context on which we 
focus our study. In the third section, we detail the derivation of our hybrid scheme which is 
based on a general conservation laws. The section 4 is devoted to the numerical results and 
comparison of our hybrid method against the WENO and ADM schemes. 


2. Biological models 

In cell biology as in physics of particles, the evolution dynamics of a group of cells or 
macro-particles in cell culture or in a bath of micro-particles plays a central role in the 
understanding and the explanation of some physical and biological behaviors. Often the 
observed quantities evolve by growth process either by nutrients uptake (the example of the 
micro-organism Daphnia which uses the nutrients to grow [15, sect. 4.3.1], [12, 16]) or by 
earnings micro-particles by polymerization (for polymers modeling [11]). 

Lot of conjectures are based on the observation of these quantities, especially on their 
evolution dynamics in long term. In modeling point of view this long term evolution dynamic 
is obtained by the analysis of the asymptotic behavior of the considered quantity. 

Following the processes taken into account in the model, the asymptotic behavior can 
either be dependent or be independent of the initial distribution of the considered group 
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of cells or parasites. Indeed, in the case where the considered population evolves only by 
growth it is proven in Lifshitz-Slyozov equations that the asymptotic behavior depends on 
the initial distribution [1, 5]. However, when aggregation process or division process is taken 
into account, the asymptotic behavior is regularized towards a quasi-universal profile as shown 
in [5] and then it is independent to the initial distribution. 

A crucial point linked to the modelling of these phenomena is the numerical simulation which 
can lead, following the used scheme, to bad conjectures on the behavior of the model. These 
bad conjectures result to some numerical artifacts caused by numerical diffusion inherent to 
some standard schemes [5]. 

In order to apply the hybrid method proposed later in this paper, we consider the following 
test case where a size-structured cells population model is taken into account and the cells 
evolve by gain or loss of micro cells. Let denote by /(i, x,£) the size density repartition of 
cells of size £ > 0, located at position x e II c M* (i = 1 , 2 , 3 ) at time t > 0 where 17 is a smooth 
bounded domain. Then the model can be written for all (t, x,£) e M + x 17 x M + as follows 


( 2 . 1 ) 


/ o o 

+ —((a(Z)c(t,x)-b)f(t,x,Q) = 0, 
f(t = 0,z,£) = /°(x,£). 


Where c(t,x ) is the concentration of micro-organisms (nutrients) and it follows a diffusion 
equation of this form: 


( 2 . 2 ) 


( d r °° 

— (c(7,x) + J o £f(t,x,£)d£) = A x c(t,x), t> 0, a; e 17 , 


—c = Vc-u = 0, on 917. 
au 


For this kind of coupling model (2.1)-(2.2), the kinetic coefficients a(£), b are interpreted as 
the rates at which cells gain or loss nutriments (monomers, micro-organisms). 

We assume the micro-cells (or monomers) to follow a diffusion equation as depicted in 
equation (2.2). We endowed this diffusion equation with a homogeneous Neumann boundary 
condition where v is the outward unit vector at point x e <917. 


The problem (2.1)-(2.2) is a variant of the very known standard Lifshitz-Slyozov sys¬ 
tem which models the evolution of a population of macro-particles immersed in a bath of 
monomers [11]. 

The analytical study of (2.1)-(2.2) concerning the existence, uniqueness and properties of 
the solution is rigorously done in [6] and the main result is based on the following hypothesis: 

Hypothesis. The kinetic coefficients a,b are required to satisfy b = 1; a is an increasing 
function with a(0) = 0 and a(+oo) = +oo; a e C°([0, oo))nC ,1 ((0, oo)) and for any £o > 0 there 
exists L a , o > 0 such that 0 < a'(£) < L a , o for £ > £o > 0. 

The initial condition satisfy c(t = 0, x) e L°°(l7); f(t = 0, x, £) € L°°(l7; L 1 ((0, oo), (l + £)d£))). 

With this previous hypothesis, the authors in [6] prove the following statement on the 
well-posedness of (2.1)-(2.2): 

Theorem 2.1. There exists a weak solution (c, /) of (2.1)-(2.2) with, for any 0 < T < oo, 
c e L°°((0, T)xl7) n L 2 (0, T; iL 1 (17)), / e L°°((0, T)xf7; L 1 ((0, oo), (l+£)d£)), c e C°([0, T]; L 2 (f7)- 
weak), / e ^([O,T]; L x (17 x (0, oo)) - weak). 
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In addition they prove thanks to the Neumann boundary condition, the following mass 
preservation relation: 


(2.3) 


d_ 

dt 


€f(t,x,€)d£dx + L c(t, x ) dx 


The fact that the space variable x acts as a parameter in the size density repartition 
/ implies that (2.1) is a transport equation and the study of its asymptotic behavior is 
numerically very challenging. 

Indeed, following the chosen model as in (2.1)—(2.2), one needs in the modelling and sim¬ 
ulations to recover the evolution dynamics of the considered population such as the time 
evolution of the total number of individuals (even in asymptotic time), the total mass of the 
population or the conservation law fulfilled by the model. For the numerical simulations of 
the evolution dynamics, a very adapted scheme is required in order to capture in exact way 
the solution of the system without artifact in order to get the right and essential properties. 
That’s the aim to introduce the following hybrid method. 


3. A HYBRID FINITE VOLUME METHOD FOR ADVECTION EQUATIONS 

3.1. Anti-dissipative method. In this section, we consider the following advection equa¬ 
tion 


(3.1) 


df_ + d(Vf) 
dt dx 


0, t > 0. 


where V(t,x ) is a given smooth velocity field. Let consider a regular mesh, with constant 
step Ax > 0: the cells are the intervals [xj_i/ 2 ,^i+i/ 2 ]) * € N with x_i/ 2 = 0, x i+ i/ 2 = (i + l)Ax, 
and x t denotes the midpoint of the cell: X{ = (i + 1/2)Ax. We denote by /” the numerical 

unknown, which is the approximation of /” = —— f 1 f(t^ n \x)dx, where = 0 < t ^ < 

Ax 

... < f( n ) < /( n+1 ) are times discretization with a possible variable time step At = t( n+1 ) - t^ n \ 
We denote by V^”jy 2 the approximations of the velocity at the cell interfaces: namely, we set 


^-!/2 = n* (n Wi/2), n € N, i € N. 


The finite volume scheme applied to (3.1) gives the following approximation 


(3.2) 


rn+1 _ nr i 

J % — J i 


Af( n ) 

Ax 


(VZl/2f? + H2-VZ ll2 fl l/2 ). 


Then the main task is how to define the range of the interface fluxes with respect to 

the stability, consistency and positivity constraints of the scheme. 

In order to describe the different constraints, let introduce the following useful notations: 


• v = 


A t^ 

Ax 


<i/ 2 = min «i)> and AP +1/2 = ma/<, 
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* * > 0 : 


1 


0 + 1/2 


rjn 

n i+ 1/2 


(ff - max(Z", ./■",)) + max(/", /” x ) 


^+1/2 

(ff ~ 1/2) + ^-1/2’ 


i+1/2 

1 


i+1/2 


(ff ~ min (ff, ff 1 )) + min (ff, /”i) 


u n V. 


- ^- 1 / 2 ) + m h/2i 


0 + 1/2 


min I B 


v: 


i+l/ 2 ’ m i-l /2 yn T ,,nm 

1 + 1/2 ^ i+ 1 / 2 , 


i— 1/2 


fi 


B 


i+1/2’ 


if m"_ |/2 > 0, 
otherwise, 


if T/ n T/ n <- D • 

11 0+1/2’ O-1/2 < u • 


h n 

0 - 1/2 


(ff - max(/", /” 0) + max(/", /” x ) 


^l^i- 1/2 

1/"|V" . | (A ? " f^i+ 1 / 2 ) + fi^i+1/2’ 


S' 


i—1/2 1 

1 


i— 1/2 


m 

0 - 1/2 


An I V-L 
0 - 1/2 


(ff ~ min (ff , /£.!)) + min(/”, 


z/ n | Vi" 

I/"IV.™ , I (/i _m i+1/2) + T 7 Z i+l/2’ 


i-1/2 

min ( 


nn 

1/2’ 




l/ 2 ’ m i+l /2 |yn | ’’’ zyW|yn 


i+ 1/21 


ff 


i- 1/21 


0 - 1 / 2 ' 


if <+ 1/2 


> 0 , 


otherwise, 


• if f^+ 1 / 2 ’ ^i- 1/2 d° not have the same sign, we set 6” +1 y 2 = &f +\/2 = ff if Vf+ 1/2 > d anc i 
^+ 1 / 2 =^ 1 / 2 = /r+i ^ ^ 1/2 < 0 . 

• /&1/2 = max «+i/2’ 6 r+i/ 2 )’ and -C1/2 = min ( M f+ 1/2,^r+1/ 2 )- 

3.1.1. Stability constraints. The stability constraints is based on the standard Courant-Friedrichs- 
Levy (CFL) condition: 


(3.3) 


0 < 


Ax 


max(|F” 1/2 |) < 1. 


From this CFL condition, we consider the case where the velocity is positive in the cell i: 
V 01/2 > 0 and V)™^ > 0, so we have 

1 


u n T0 


- 1 > 0 


i+1/2 


and in first hand we write 

Then we deduce 
(3.4) 


>«y n , 

i+1/2 


!)(/"- m in(/f,/£i)) > 0. 


^ v i+1/2 


(/I - min(/”,/”,)) + min(/", /",) > /” 
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In second hand we have 

I 

then we deduce 
(3.5) 


' v n V n 

i+l/2 


i)(/r 


max 


(/r,/r_i))<o, 


h , n V n , 
u v i+l/2 


(/” - m ax(/'", /",)) + max(/”, /"j) < ff 


In last hand, assuming that /" > 0 Vz (the positivity constraint is discussed later) we deduce 
the obvious relation 

V n f r ' 

min(/f,/f_ 1 )-t^ + 


/T 


p u n V n v n V n 

V i+l/2 V V i+ 1/2 ^ *+1/2 


(3.6) 

' ^ fl¬ 
owing to (3.4)-(3.6) and the previous notations, one deduce this first non empty interval 
where the adequate flux will be chosen 

(3-7) f^[bl i/ 2 ,#ki/ 2 ]*0. 

For the case where the velocity is locally negative mean VP/y 2 < 0 and F™-w 2 < 0, we 
perform the same reasoning and obtain the non empty interval /” e [^ 1 / 2 ’^^- 1 / 2 ] * 


3.1.2. Consistency constraints. For the consistency of the scheme, we write the very definition 
which is that the numerical flux between cell i and cell i + 1 belongs necessary to the interval 
defined by the numerical solutions /” and /”. j. So the consistency constraint is written as 
follows 


(3.8) 


m ti/2 ^ fi+1/2 * < + i/ 2 , Vi in the grid. 


3.1.3. Positivity constraints. From a positive initial solution, we need to impose a condition 
on the numerical fluxes in order to ensure the positivity of the numerical approximation of 
the solution given by (3.2). For this, we are looking for conditions so that 


/r 1 * 0 


V n W+l/ 2 fi+l /2 - K 1 / 2 &I/ 2 ) * fi 


rn - rn 

J i+l/2 ~ h- 


V n 
V i- 1/2 

1/2 yn 

i+l/2 


rn 

J i 


u n V n , 

i+l/2 


so, obtaining the later inequality is done by imposing the flux to satisfy 


(3.9) 


yn 

+ 1+1/2 1-1/ 2 yn 

i+l/2 


Vi 


jyTl'yn 

i+l/2 


because by the consistency constraints we know that vn!f_ Y j 2 < /” 


Proposition 3.1. Assuming the CFL condition (3.3) be satisfied. Then for any i the interval 
[^+ 1 / 2 , ^ 1 / 2 ] is non empty- B y choosing the fluxes /” 1/2 e [/r™ +1/2 , ^+i/ 2 ] f or an V then 
the following assertions hold: 

1) The scheme (3.2) is consistent with (3.1). 

2) The scheme (3.2) remains positive for positive initial solution: mean if ff > 0 for any 
i then ff +1 > 0 too. 

3) The scheme satisfies: if Vff > 0 then mHf_y 2 < ff +1 < MP_y 2 , while if Vf 1 < 0 then 


m i+ 1/2 


< f n+1 < M n 
-n - m i+ 1 / 2 - 


Proof. The proof of the proposition is essentially based on the gathering of the results (3.4)- 
(3.6) and (3.8)-(3.9). For more details on the proof, one can refer to [5]. 
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(3.10) 


T+l/2> 

* ft 1 * 


'/n 

u+1/2 ’ 


Remark 3.2. In the previous reasoning, we exclude the case where the velocities on the 
interfaces cell are of different sign. So in the case where Vfl j y 2 < 0 and VIf 1 j 2 > 0 then it is 

obvious that the positivity of the numerical solution ff +1 is obtained without any time step 
condition. Nevertheless, the case where Vfl j , 2 > 0 and Vflij 2 < 0 mean the possible empty 
of the cell from the two sides. So we choose the upwind fluxes and the numerical solution 
becomes constant in the cell i. 

Remark 3.3. In this above presentation, the numerical flux is designed for general conser¬ 
vation laws, while the one in [5] is more suitable for the transport equations. 

For the anti-dissipative strategy, we define the flux ffl +1 / 2 by solving the following mini¬ 
mization problem: 

To minimize | / t " 1/2 - f? +1 \ 
under the constraint /” 1/2 e [^ +1/2 , ^+i/ 2 ] 
which solution is given for instance in the case V^” 1 , 2 > 0, V- n > 0 and > 0 by 

f n - n n if f n < n r 

h+ 1/2 “ ^i+l/2 11 H +1 - h'i 

fl 1/2 = K i ^ A&i, 

/r + i /2 =^ n + i /2 ^ /r + 1 * 

This kind of anti-dissipative method is very suitable for discontinuous initial solution, 
which has been shown in [4], However, it is not suitable for smooth solution. Indeed, it turns 
the very smooth solution into a series of step functions with respect to time evolution [4, 5]. 
The objective of remain part is to find an alternative method such that it keeps the shock 
near discontinuities while has high accuracy in smooth regions. To the end, we propose the 
following hybrid method. 

3.2. A hybrid method. We denote /A 1/2 the flux computed by the anti-dissipative method 
and f™ 1/2 the flux computed by a high accurate method, for instance the fifth order method [8]. 
So our desirable flux f^h/ 2 by the hybrid method will just be a convex combination of , 2 

and 1/2 > ie - 

(3-11) fl 

where = 1, 1/2 - Moreover, it is desirable: 

* “ti /2 = 0(1) and w,W /2 = o(l), near discontinuities, 

• = 0(Ax ai ) and ^Jh/ 2 = 0(1), in smooth regions, 

where M is a large enough number that causes f^±/ 2 1° be the dominant term in smooth 
regions. 

According to these considerations, we then need to identify the smoothness of solution. 
We consider a similar smoothness measurement, which was first proposed in [3] 

(3.12) „ 

where D is a scaling value given as 

(3.13) 


_,, A f A ,,,w ,w 

“ U i+l/2h+l/2 + W j+l/2/j+l/2 


e,: = 


D 


D = max fl - min fl 


Figure 1 shows the interpolated value, i.e. fi is obtained by the fourth-order interpolation. 
Note that the point fl itself, as shown in Figure 1, is not included in the interpolation. 


fl - ~(-fl~2 + 4/j-i + 4/j + i - /j +2 ). 
0 


(3.14) 
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Figure 1. Interpolation error between the node value fi and interpolated value fi. 


Then we choose the smoothness measurement at the cell interface by an upwind way 

e i-> if V i+ 1/2 > 0, 
ej+i, else. 


(3.15) 


°i+ 1/2 


Finally, the weights have forms 

(3.16) w£ 1/2 = 1 - exp(-e- +1/2 /c x ), ug 1/2 = exp (~e 2 i+1/2 /c x ), 

where 

<-) -(f)"- 

The parameter a is a constant, and a discussion of its value will be given in the section 4.1.1. 
We now present the order analysis of the weights. In smooth regions, using a Taylor expansion 
analysis, we have: 

e i+1/2 = 0( Ax 4 ). 

Now using the Taylor expansion of the exponential function and assuming a < 1, one obtains: 


uj; 


i+i /2 = 1 - e M-e i+1/2 /c x ) * e i+1/2 /c x = 0(Ax 8 a ), 


and 


W i+l/2 _ ex p( e i+l/2l Cx ) ~ 1 e i+l/2/ C z _ 0(1). 


°i+l/2/ Cx ) w f e i+ 

Therefore, we use the fifth order WENO flux in smooth regions, i.e. 

rH „ f w 
1 i+l/2 Ji+ 1/2■ 

For a discontinuity, since the size of the discontinuity does not change as Ax ->■ 0, one can 
conclude 

e i+l/2 = 0(1). 

Therefore, applying the same analysis used above leads to 

Al/2 = 0(1). 

Thus, at discontinuities the anti-dissipative flux is activated. 


Remark 3.4. In the original method [3], the parameter D depended on the index i and it 
was a convex combination of the global and local scales, i.e. 

Di = C S Sg + (1 — Cg') Si, 

where S g = max*./) - minj fi, i e {1,..., w x }, S t = max* f - min* /*, i e {i - 2,..., i + 2}. The 
parameter c s was chosen as 0.1 or 0.01 in [3]. This choice can highlight the small jumps in 
the solution. However, in the paper we focus on the major jumps in the solution. Therefore, 
the global scale S g seems more appropriate. 
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Remark 3.5. Again, in the original method [3], the parameter c x is a constant, which varies 
from problem to problem. At contrast, we propose the parameter c x depending on the scaled 
mesh step, i.e. Ax/L x and the exponent a. Moreover, according to the numerical experience 
in section f. 1 . 1 , the parameter a can be a fixed number. 


4. Numerical results 

In this section, we will present several numerical results to depict the behaviors of our 
hybrid method and its applications in population dynamics. 

4.1. The classical numerical tests. Here, we first discuss the parameter a in (3.17), which 
is the key point in our hybrid method. Then we perform the classical tests with the ID free 
transport equation and the 2D rotation equation to verify the convergence of our hybrid 
method for the regular and irregular initial data. In the sequel, the third order Runge-Kutta 
method is used for time discretization. 


4.1.1. A discussion of the parameter a in (3.17). To determine the parameter a, we use the 
free transport equation 


(4.1) W + i£ = 0, ^°> 

with a very oscillatory initial condition, given by [8] , 

fi(x) = jj[G(x, z- 5) + G{x, z-5) + 4 G(x,z)], 
f2(x) = 1 , 


(4.2) 


/( 0 ) x) = - 


h{x) = l-|10(x- 0.1)1, 

U(x) = \[F(x,z- 8 ) + F(x, z-5) + 4 F(x,z)], 


0 , 


if - 0.8 < x < -0.6, 
if - 0.4 < x < -0.2, 
if 0 < x < 0.2, 
if 0.4 < x < 0.6, 
otherwise. 


where G{x,z) = exp (~/3(x - z) 2 ), F{x,a) = {max((l - a 2 {x - a) 2 ) 1//2 ,0)} with a = 0.5, 
z = -0.7, 5 = 0.005, a = 10 and /3 = (log2)/36d 2 . Moreover, the periodic boundary conditions 
is considered. 

The numerical results is presented in Figure 2. We first see that the WENO method is 
well adapted for the smooth regions of solution, while it becomes significantly diffuse near 
the steep regions. For the refined mesh, we can always observed this diffusion near the step 
function clearly. The anti-diffusive method has perfect performance for the step function, 
however it destroys the very smooth function. Indeed, it turns the smooth function into a 
sawtooth profile, which can not be regularized by refining mesh. 

Now we expect that the hybrid method can combine the advantages of both the WENO 
method and the anti-diffusive method. However, it remains to choose an appropriate param¬ 
eter a in (3.17). In Figure 3, we present the numerical solutions of the hybrid methods with 
the different values of the parameter a , i.e. a - 0.65, 0.75, 0.85. We first observe that the 
hybrid method approximates well the functions f \ and fj, of (4.2) for all the values of a. How¬ 
ever, there are significantly differences for the function /2 and f±. Indeed, with a = 0.65, the 
hybrid method can not preserve the steep profile of the step function /(>; while with a = 0.85, 
the sawtooth profile appears again in the function f 4 . By refining the mesh, we see that 
the approximations of the functions f 2 and f± are improved, but they are still not perfect. 
Finally, with a = 0.75, we are capable to capture well all configurations of the solution (4.2). 
Therefore, we will use the value of the parameter a = 0.75 in the sequel of this paper. 
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Figure 2. Numerical solutions of the WENO method and the anti-diffusive 
method for ID transport equation : Plot solutions of the linear equation (4.1) 
with initial data (4.2). n x = 200 in the first row and n x - 400 in the second 
row, while At is chosen such that CFL- 0.2. The final time is T en d = 8. 



Figure 3. Comparison of different values of the parameter a in (3.17) of 
the hybrid method for ID transport equation : Plot solutions of the linear 
equation (4.1) with initial data (4.2). n x = 200 in the first row and n x = 400 
in the second row, while At is chosen such that CFL- 0.2. The final time is 
Tend = 8 . 
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4.1.2. Convergence rate for the ID transport equation. Let us first consider a smooth solution, 
where the initial condition is chosen as 

/(0, x) = sin (irx) , xe [-1,1]. 

The numerical error for different methods is presented in Table 1. We first notice that 
the anti-diffusive method has only first order of convergence rate. At contrast, the WENO 
method and the hybrid method have both the fifth order convergence rate. Thus they are 
much more precise than the anti-diffusive method. Moreover, the hybrid method does not 
disturb at all the precision of numerical results in smooth solution case. 


Tlx 

200 

400 

800 

|| • ||i 

r 

|| • ||i 

r 

|| • 111 

r 

Anti-diffusive method 

4.31e-2 

1.07 

2.06e-2 

1.07 

1.06e-2 

0.96 

WENO scheme 

1.14e-7 

5.00 

3.57e-9 

5.00 

1.12e-10 

4.99 

Hybrid method 

1.14e-7 

5.00 

3.57e-9 

5.00 

1.13e-10 

4.99 


Table 1. ID transport equation : Error in L\-norm and order of convergence 
r. The final time is T en d = 8. 


We next consider a step function as follows 

[ 1, for - 1 < x < 0, 

(4.3) /(0,s) = j 

[ 0, otherwise. 

Comparisons of the methods are now summarized in Table 2. We observe that the anti- 
diffusive method has the first order convergence rate for irregular solution, while the WENO 
method can not achieve the first order convergence rate. The hybrid method is less precise 
than the anti-diffusive method, however it is more precise of one degree than the WENO 
method. It is interesting to notice that the hybrid method has a convergence rate even more 
than one for irregular solution. At last, we see that the total variation for the hybrid method 
is comparable with the WENO method, which means that the hybrid method is also an 
essentially non-oscillatory method. 


n x 

200 

400 

800 

|| • ||i 

r 

|| • ||i 

r 

|| • ||i 

r 

Anti-diffusive method 

1.48e-3 

1.00 

7.39e-4 

1.00 

3.69e-4 

1.00 

WENO scheme 

4.50e-2 

0.83 

2.53e-2 

0.83 

1.43e-2 

0.82 

Hybrid method 

7.93e-3 

3.32 

2.36e-3 

1.75 

9.24e-4 

1.35 


(a) Error between exact solution and approximated solution 


n x 

200 

400 

800 

Anti-diffusive method 

-1.33e-15 

-1.55e-15 

-1.78e-15 

WENO scheme 

9.71e-7 

9.52e-7 

7.75e-7 

Hybrid method 

8.51e-7 

5.28e-6 

1.16e-6 


(b) Error of total variation 


Table 2. ID transport equation : Comparison of different methods for the 
linear equation (4.1) with initial data (4.3). (a) Error in L\ norm and r is the 
order of accuracy (b) Error on the total variation. The final time is T en <i = 8. 
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Finally, the hybrid method does not significantly increase the computational cost. Indeed, 
in average the WENO scheme and the ADM method take 74 operations and 27 operations 
respectively. While the hybrid method consists of these two methods (103 operations) and 
in addition the computation of smoothness indicator (27 operations). So in total the hy¬ 
brid method takes approximately 130 operations, about 1.76 times of the WENO scheme or 
4.81 times of the ADM method. However, considering its important features, we think the 
additional computational cost of our hybrid method is acceptable. 


4.1.3. 2D rotation equation. We can directly extend the hybrid method to two dimensional 
cases. In fact, the numerical flux can be computed dimension by dimension. For example, in 
the x direction, the numerical flux can be computed by 

(A A) fH _ A j-A w ,w 

l 4 ' 4 ! H+l/2,j - U i+l/2,jh+l/2,j + U i+l/2,jh+l/2,j 

where oj A . . . = 1, ur 4 , /Q ., /0 . > 0. Then we choose the smoothness measurement 

1+1/2,j 1+1/2,j ’ 1+1/2,j ’ t+l/2,j 

at the cell interface by an upwind way 


(4.5) 


e i+l/2,j 


hr 


if k)+i/2j ^ 0) 


e 


X 


else, 


where the smoothness measurement e® ■ is the same as in (3.12)-(3.14) for j fixed. The weights 
have forms 


(4.6) 
where 

(4.7) 


U tl/ 2 ,j = 1 - ex p( -e i+ 1/2 j/ c z), ^T+l/ 2,3 = exp (-e 2 i+ii 2 j/c x ), 



Similarly, we can compute the numerical flux /^ +1 / 2 i n V direction. 

To test the efficiency of the hybrid method, we use the famous Zalesak’ disk test. The 
governed equation is the 2D rotation equation 


(4.8) 


df 

dt 


+ y 


df 

dx 


The initial solution is [9] 
(4.9) 


dj_ 

X dy 


f(0,x) = ■ 


l-^\/x 2 + (y + 0.5) 2 , 

1+cos 7T y- (x+0.5) 2 +y 2 
4 ’ 


0 , 


= 0, x x y € [~1,1] x [—1,1], t> 0. 


if ^x 2 + (y- 0.5) 2 < r 0 && (|x| > 0.025 || y > 0.75), 
if \/x 2 + (y + 0.5) 2 < ro, 

if yj(x + 0.5) 2 + y 2 < ro, 
otherwise, 


where the radius ro = 0.3. 

We illustrate the numerical results in Figure 4. We first notice that the initial condition 
consists of the Zalesak’ disk, the conical body and the peak of the hump. Then after one 
period, the WENO method preserves well the the conical body and the peak of the hump, 
however a clear diffusion appears in the Zalesak’ disk. The anti-diffusive method keeps well 
the shape of the Zalesak’ disk, but it destroys completely the others two objects. Finally, 
we observe that the hybrid method performs well for all these three objects. Therefore, the 
hybrid method is suitable for both the smooth solution and the irregular solution for the 
advection equations. 
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Contour 
Var: density 


1 

: 


0.947368 

0.894737 

0.842105 

0.789474 

0.736842 

0.684211 

0.631579 

0.578947 

0.526316 

0.473684 

0.421053 

0.368421 

0.315789 

0.263158 

0.210526 

0.157895 

0.105263 

0.0526316 


Max: 1.000 
Min: -0.04013 



(a) Reference solution 


(b)WENO method 


Var: density 


-0.947368 -i 

■ 

-0.894737 


-0.842105 - 


0.789474 


0.736842 - 


-0.684211 


-0.631579 - 

■ 

-0.578947 


-0.526316 - 

■ 

0.473684 ; 


0.421053 ; 

| 

-0.368421 

1 

-0.315789 . 


-0.263158 


-0.210526 J 


-0.157895 

■ 

0.105263 J 


0.0526316 


-0 



Max: 1.000 
Min: 0.000 


Contour 
Var: density 

1 0.947368 
0.894737 
0.842105 
0.789474 
-0.736842 

1 0.684211 
0.631579 
0.578947 
0.526316 
(-0.473684 
-0.421053 

■ 0.368421 
0.315789 
-0.263158 
(0.210526 
-0.157895 

f 0.105263 
0.0526316 
■ 

Max: 1.000 
Min: -0.0001130 



(c) Anti-diffusive method (d) Hybrid method 

Figure 4. 2D rotation equation : Plot solutions of the linear equation (4.8) 
with initial data (4.9). Mesh size is n x x n y = 400 x 400, while At is chosen 
such that CFL- 0.2. The final time is T en d = 27t. 


4.2. Biological numerical tests. In this subsection, we will focus on numerical simulations 
for the polymerization/depolymerization type model. The objective here is to highlight the 
good performance of our hybrid method on the long term asymptotic behavior of the solution. 


4.2.1. Space-homogeneous polymerization/depolymerization type model. Here we are inter¬ 
ested on the numerical behavior of the standard Lifshitz-Slyozov equations [11] that we refer 
as homogeneous in space in comparison to the system (2.1)-(2.2). 

This standard Lifshitz-Slyozov equations can be interpreted in the point of view of popu¬ 
lation dynamics as a model describing a population of cells evolving only by nutrients uptake 
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Figure 5. Polymerization/depolymerization test in homogeneous space: the 
regular and irregular initial data corresponding to (4.11) and (4.12). 


(without birth, death and division) where the nutrients are characterized by their concentra¬ 
tion c(i) which fulfills a mass preservation equation as follows 


(4.10) 


J)/(«) + |l« 1/3 c(t) -l)/(t,0) = 0, 

J f' oo 

£f(t,Od£ = p > 

0 


where p is a constant and measures the total initial mass, £ the cell-size and / the size density 
repartition. 

Despite its simplistic appearance, this model (4.10) is very intriguing when one is interested 
on the time asymptotic behavior. An interesting discussion is made in [5] and the authors 
highlight specifically the importance of using an anti dissipative numerical scheme in order to 
avoid numerical artifact (diffusion) and then capture the exact asymptotic profile. In order 
to test our scheme defined in section 3 and show its accuracy, we compare the numerical 
results with those obtained either by a WENO scheme (see [1]) or by the anti-dissipative 
scheme (see [5]). Two types of initial distribution functions will be considered: the regular 
one, which represents the size-density of cells in a normal distribution as in (4.11), 


(4.11) 


f°(x) = 0.1 exp(-0.1(£ - 20) 2 ). 


While the irregular one represents that the size-density of cells is concentrated around a spot 
as in (4.12), 


(4.12) 


no 


1, if 10 <£<30, 

0, else. 


These two types of size-density are illustrated in Figure 5. 

We first consider the numerical results with the regular data (4.11) (see Figure 6). Our 
hybrid scheme has the same results as the classical WENO scheme. More precisely, a very 
smooth profile is formed, which moves towards cells of large size. This observation has a 
good agreement with the physical behavior of transport equation which does not modify the 
general shape of the initial distribution function. In biological point of view, the observation 
show the growth of big size cells at the expense of smaller ones. This competition between 
big size cells and smaller size ones is due to the fact that at any time there exists a critical 
size £ C rit = such that the velocity vanish; then cells of size £ > £ cr jf grow while cells of 
size £ < £ cr it shrink. This competition phenomenon is well known under the name “Ostwald 
ripening” [13, 14, 10, 17]. At contrast, the anti-dissipative method forms a very sharp front, 
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which is caused by the Downwind flux and it can not be significantly improved by mesh 
refinement. Even in the zone view for small size particles (see the right column of Figure 6), 
the anti-dissipative method generates “stairs” looking like oscillations. 

Then, with the irregular data (4.12), as pictured at the figure 7, we first notice that the 
irregular data is regularized by the numerical diffusion by using the WENO scheme. Indeed, 
I. M. Lifshitz and V. V. Slyozov in their original paper [11] conjecture that the asymptotic 
solution is in the same form despite the initial data (compare with the left column of Figure 6). 
However, as pointed out in [5], the irregular data does have an influence for long term solution. 
Our scheme as the anti-dissipative one preserves well the front propagation at long time (no 
numerical diffusion) what is physically very important when dealing with transport equation 
and can be interpreted in biology as the fact that very stiff localized population of cells 
must remain stiff localized if they are subject to only transport process. Moreover, our 
scheme corrects well the “stairs” looking like oscillations in the numerical solution of the 
anti-dissipative method (see the right column of Figure 7). 


4.2.2. Non-space-homogeneous polymerization/depolymerization type model. Here we are in¬ 
terested on the numerical behavior of the non-homogeneous system (2.1)-(2.2) with the choice 
a(£) = £ 3 / 3 , b - 1, that we rewrite as follows 


(4.13) 


/ Q O 

—f(t,x,0 + — ((C 1 / 3 c(t,x) - 1 )f(t,x,0) = 0, t > 0 , x € n, £ > 0 , 
d t [c{t,x) + £f(t,x,£)d£j = A x c(t,x), t> O.nefl. 


The diffusion equation is endowed with homogeneous Neumann boundary condition 


d v c = Vc • v = 0, on <9D. 

Here we choose fl c M means ID in space-variable. The system (4.13) describes the immersion 
of a population of cells in a culture of micro-organisms (nutrients) that are subjected to 
diffusion equation while de size density repartition of cells is parametrized by the space 
position x. 

The numerical simulations are performed in the slab x e [0,60]. The size variable is 
truncated to £ € [0,100]. The initial concentration is defined by 

(4.14) c(t = 0,x) = 0.5I x6 [ 20) 4o] • 

As we did in the previous test, two types of size-density are considered: a regular one 

(4.15) f(t = 0, x, 0 = 0.01 exp(-0.2(£ - 30) 2 )I xe[20 , 40] . 
and an irregular one 

(4-16) f(t = 0,X,£) = 0.01I xe [20,40] x %e[30,35]- 

The initial data is presented in Figure 8. 

Let us briefly describe the algorithm for the system (4.13). The third order Runge-Kutta 
method is used as the time discretization for the advection equation of the size-density of 
macro-particles. In the reconstruction of numerical flux, both the WENO scheme and our 
hybrid scheme are applied for the purpose of comparison. Then for the diffusion equation of 
monomers, the crank-Nicolson method is used for time discretization. A classical second-order 
finite volume method is applied for the lapacian operator, and the integral is approximated 
by the classical Simpson’s rule. 

Firstly, in the case with the regular data (4.15), there are almost no difference from the 
size-density distribution between the two numerical methods (see Figure 9). The initial 
distribution is in a plaque form at t = 0, then it generates a moon shape caused by the 
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Figure 6. Polymerization/depolymerization test in homogeneous space: Plot 
solutions of the equation (4.10) with the regular initial data (4.11). Mesh size 
is n x = 800, At = 1/8, CFL » 0.12, p = 41. 
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(a) t = 1250 




(b) Zone, t = 1250 

x 10" 4 



(c) t = 1500 


(d) Zone, t = 1500 




Figure 7. Polymerization/depolymerization test in homogeneous space: Plot 
solutions of the equation (4.10) with the irregular initial data (4.12). Mesh 
size is n x = 800, At = 1/8, CFL & 0.12, p = 41. 
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X-Axis 
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0.004996 

0.002498 
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(a) Regular initial size distribution function (b) Irregular initial size distribution function 



(c) Initial concentration 


Figure 8. Polymerization/depolymerization test in non-homogeneous space: 
the initial data corresponding to (4.14)-(4.16). 


diffusion of the concentration of monomers. In evolution in time, we see that the solution is 
very smooth in size direction. 

Secondly, we consider the irregular data (4.16). The micro-organisms (nutrients) concen¬ 
tration c(t,x ) and the mass of the cells / 0 °° £/(£,£, £)d£ are presented in Figure 10. From 
these two quantities, we see there are no difference at different times. Moreover, the total 
mass for both methods preserves well in whole time evolution (see Figure 11). However, we 
observe clearly that numerical dissipation appears with the WENO scheme, while our hybrid 
method preserves a sharp profile (see Figure 12). 

In general, having a scheme able to reconstruct accurately the profile of the solution is of 
paramount importance in population dynamics; indeed this profile is often used for parameter 
estimation and comparison with experimental measurements for prediction purposes such as 
prediction of the evolution of bacterial populations followed by genetic algorithms or the 
assessment of the average time for the balance of a product (micro-organisms or monomers 
in our case) in the considered mixture. That’s the case in this paper because following the 
initial distribution, the expected solution does not have the same profile and thus does not 
lead to same predictions. 
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Figure 9. Polymerization/depolynrerization test in non-homogeneous space: 
Plot the size distribution function of the equation (4.13) with the regular initial 
data (4.14), (4.15) at different time. Mesh size is n x xn^ = 100x800, A t = 0.1, 
CFL Ri 0.13. 
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X axis 

(a) WENO scheme 


X axis 

(b) Hybrid Method 




(c) WENO scheme 


(d) Hybrid Method 


Figure 10. Polymerization/depolymerization test in non-homogeneous 
space: Evolution of the monomers concentration c(t,x ) (first row) and the 
mass of the marco-particles / 0 °° £f(t,x,£)d£ (second row) corresponding to 
the equation (4.13) with the irregular initial data (4.14), (4.16). Mesh size is 
n x x = 100 x 800, At = 0.1, CFL & 0.13. 



Figure 11. Polymerization/depolymerization test in homogeneous space: 
Mass conservation property corresponding to (2.3). Mesh size is n r x nc = 
100 x 800, At = 0.1, CFL * 0.13. 
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Figure 12. Polymerization/depolynrerization test in non-homogeneous 
space: Plot the size distribution function of the equation (4.13) with the irregu¬ 
lar initial data (4.14), (4.16) at different times. Mesh size is n x xn^ = 100x800, 
At = 0.1, CFL 0.13. 
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5. Conclusion and perspective 

In this paper we have proposed a hybrid finite volumes scheme based on the flux convex 
combination between the Anti-Dissipative Method (ADM) [4, 5] and the WENO order 5 
method [8]. The obtained numerical results show a good accuracy in reconstructing the 
solution of transport type equation even in case of discontinuous initial data. Indeed the 
simulations in figure 6 show results that are as good as the ones in WENO scheme and better 
than ADM scheme which fail for regular initial distribution. In reverse, for irregular initial 
distribution, the hybrid scheme show better numerical results than de ADM scheme in the 
sense that it advects very well the solution without “stairs” and it shows also better results 
than WENO scheme which develop numerical diffusion for irregular distribution as depicted 
in figure 7. So, when the WENO order 5 scheme fails because of numerical diffusion artifact, 
the hybrid method remains anti-dissipative and when ADM scheme develops “stairs” like 
oscillations, the hybrid scheme corrects them. This property is very suitable for long term 
asymptotic behavior of the solution of population dynamics, as presented in the numerical 
simulations for the polymerization/depolymerization type models. 

Finally, lets point out that the flux construction of our hybrid scheme depends on an empir¬ 
ical parameter a= | which is related to the weights in the convex combination. Nevertheless 
we can’t explain yet why exactly a = | but we are working in a preprint paper which is 
focused on the theoretical and numerical explanation of this parameter a. 
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